/// Alvaro Cuno, Mar2008
/// MatrixT definition file

#ifndef _MatrixT_H__
#define _MatrixT_H__

#include <iostream>
#include <string>
#include <cassert>

using namespace std;

/// Defines a MatrixT
template <class Number>
class MatrixT {
public:
   /// Default constructor
   MatrixT(): r(0), c(0), p(0), m(r), n(c) { }
   /// Copy constructor
   MatrixT(const MatrixT &M): r(M.r), c(M.c), m(r), n(c) {
      p = new Number *[r];
      p[0] = new Number[r*c];
      for (int i=1; i<r; ++i) p[i] = p[i-1] + c;
      *this = M;
   }
   /// constructor
   MatrixT(int m, int n): r(m), c(n), m(r), n(c) {
      p = new Number *[r];
      p[0] = new Number[r*c];
      for (int i=1; i<r; ++i) p[i] = p[i-1] + c;
   }
   /// destructor
   ~MatrixT() {
      if (p != 0) {
         delete [](p[0]);
         delete []p;
      }
   }
   /// Assignment operator
   const MatrixT& operator=(const MatrixT &M) {
      if (this == &M) return *this;
      if (r != M.r || c != M.c) {
         if (p != 0) {
            delete [](p[0]);
            delete []p;
         }
         r=M.r; c=M.c;
         p = new Number *[r];
         p[0] = new Number[r*c];
      }
      for (int i=1; i<r; ++i) p[i] = p[i-1] + c;
      for (int i=0; i<r; ++i) 
         for (int j=0; j<c; ++j) 
            p[i][j] = M[i][j];
      return *this;
   }
   
   Number* operator[](const int i) {
      return p[i];
   }
   const Number* operator[](const int i) const {
      return p[i];
   }
   Number& operator()(const int i, const int j) {
      return p[i][j];
   }
   const Number& operator()(const int i, const int j) const {
      return p[i][j];
   }
   
   /// MatrixT + MatrixT operator
   MatrixT operator+(const MatrixT &B) const {
      MatrixT C(*this);
      assert(C.c == B.c and C.r == B.r);
      for (int i=0; i<r; ++i) 
         for (int j=0; j<c; ++j) 
            C[i][j] += B[i][j];
      return C;
   }
   /// MatrixT - MatrixT operator
   MatrixT operator-(const MatrixT &B) {
      MatrixT C(*this);
      assert(C.c == B.c and C.r == B.r);
      for (int i=0; i<r; ++i) 
         for (int j=0; j<c; ++j) 
            C[i][j] -= B[i][j];
      return C;
   }
   /// MatrixT * MatrixT operator
   MatrixT operator*(const MatrixT &B) const {
      const MatrixT &A = *this;
      assert(A.c == B.r);
      MatrixT C(A.r, B.c);
      C.zero();
      for (int rr=0;rr<A.r; ++rr) 
         for (int cc=0;cc<B.c; ++cc) 
            for (int i=0;i<A.c; ++i) 
               C[rr][cc]+=A[rr][i]*B[i][cc];
      return C;
   }
   // Sets all the entries to zero value
   void zero() {
      for (int i=0; i<r; ++i) 
         for (int j=0; j<c; ++j) 
            p[i][j] = 0.0;
   }
   // Sets an identity MatrixT
   void identity() {
      for (int i=0; i<r; ++i) 
         for (int j=0; j<c; ++j) 
            p[i][j] = (i==j)? 1.0:0.0;
   }
   
   int nrows() const { return r;}
   int ncols() const { return c;}
   
   /// output operator
   template <class T>
   friend ostream &operator<<(ostream &, const MatrixT<T> &);
   
public:
   int      &m; // an alias for the number of rows
   int      &n; // an alias for the number of columns
private:
   Number **p; // pointer to the entries of the MatrixT
   int      r; // number of rows
   int      c; // number of columns
};

/// output operator
template <class T>
ostream &operator<<(ostream &out, const MatrixT<T> &A) {
   for (int i=0; i<A.r; ++i) {
      for (int j=0; j<A.c; ++j) 
         out<<A[i][j]<<" ";
      out<<endl;
   }
   return out;
}

#endif
